Hydrodynamic theory for dissipative hard spheres with multi-particle interactions 
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Extension to kinetic theory and hydrodynamic models are proposed that account for the existence 
of multi-particle contacts. In the presence of multi-particle contacts (involving elastic, reversible, 
potential contact energy), dissipation of the translational (kinetic) energy is reduced and a class 
of different models lead to deviations from the classical inelastic hard sphere (IHS) homogeneous 
cooling state (HCS), as examined here. The theoretical results are found to be in perfect agreement 
with the numerical simulations. 
PACS: 45.70, 47.50+d, 51.10.+y, 47.11. +j 



o 



CO 

a 

I 

C 

O 

o 



> 



(N 
O 

a 

c 

o 
o 



X 



Kinetic theory and the related hydrodynamic models 
are helpful tools for the modeling and understanding of 
transport processes in classical, elastic gases for low and 
moderate densities (J, ^ . The hard sphere (HS) model is 
the corresponding approach to be implemented as a nu- 
merical model || Q. A successful theoretical approach 
in the spirit of Boltzmann or Chapman and Enskog jj], ||] 
requires the basic assumptions: (i) The collisions are in- 
stantaneous and (ii) subsequent collisions are uncorre- 
cted ("molecular chaos"). Conditions (i) and (ii) lead 
to the Boltzmann equation, and in the equilibrium state, 
the velocity distribution is (iii) a Maxwellian. 

When dissipation is added, one has the inelastic hard 
sphere (IHS) model, where the coefficient of restitution 
r quantifies dissipation, elastic systems have r = 1, and 
1 — r 2 > determines the amount of energy lost in a two- 
particle collision in the center of mass system. The range 
of applicability of the theory for the IHS was addressed 
in several papers ||, @, §]• However, this is far from 
the scope of this study, since it involves e.g. Sonine poly- 
nomial expansions and "viscoelastic" material behavior 
f|, so that we just assume (i-iii) valid for the sake of 
simplicity. 

In this study we restrict ourselves to the homogeneous 
cooling state (HCS) and focus on a mean-field hydrody- 
namic approach || @, [§, 0, 0, neglecting spatial 
structures like clusters or shear modes. This idealization 
is reasonable for either weak dissipation, low density or 
small system size, however, we do not discuss the range 
of applicability here. The qualitative prediction for the 
long-time decay of energy with time by Haff |lj] is con- 
firmed and it was shown that the distribution function is 
isotropic in velocity space and it is close to a Maxwellian 
as long as the system is homogeneous || [u| |i~4| . 

In the IHS collisions are always instantaneous, see con- 
dition (i), due to the rigid interaction potential. On the 
first glance, this makes the model (and kinetic theory 
too) inadequate for the description of real materials for 
which the interaction potential may be steep, but is never 
perfectly rigid, see Fig. [l] During the contact of two real 
particles, kinetic energy is stored in elastic (reversible) 



potential energy that, in the static limit, can be recov- 
ered after very long times. The conclusion is thus that 
one has a fraction of elastic energy, which cannot be dis- 
sipated, in the real system. This fraction is missing in all 
idealized models HS, IHS, and also in the kinetic theory, 
and has to be defined. Thus we will propose and exam- 
ine possible ways to cure this problem of the hard sphere 
model, but keeping kinetic theory still applicable. 
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FIG. 1: Schematic plot of the trajectories of two soft (left) 
and two hard (right) particles against time. The beginning 
and the ending of the interaction are marked by dashed and 
dotted vertical lines, respectively, and the time t c during that 
kinetic energy is stored as elastic energy in the contact - so 
that dissipation is affected - is marked as shaded region. 

The first step is to define or identify possible multi- 
particle contacts. In a real system (or in a soft-particle 
model) one just counts the number of contacts a particle 
has. Within the extented IHS model, a particle remem- 
bers its last contact and every contact occuring within a 
time t c after that, see the shaded area in Fig. [l], is defined 
a multi-particle contact. 

In low density systems, where the mean free flight time 
is much larger than the contact duration, multi-particle 
contacts are rare. However, assumption (i) can also be 
valid in high density situations where the free path is 
much smaller than the particle diameter: This is possi- 
ble in the case of an extremely short contact duration, 
when collisions remain practically instantaneous. Thus 
we conclude that the free path, i.e. the density, is not an 
appropriate measure for the occurence of multi-particle 
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contacts. We rather define, as a more objective criterion, 
the ratio r c = t c /tE between the contact duration t c and 
the typical time between collisions tg as obtained by the 
Enskog theory R ^ [l^] . Small and large r c values corre- 
spond to pair- and multi-particle collisions, respectively, 
in our framework |D| . 

Up to now, the elastic HS model is not changed at all 
concerning particle trajectories or whatsoever. The only 
modification is, that every particle that had a collision a 
short time ago keeps this event in memory for the contact 
duration t c and every new contact occuring within this 
time is now defined as an elastic contact. This allows to 
split the total energy in the system into a kinetic and an 
elastic (potential) contact energy jl6| , as can be done in a 
real system. If a part of the kinetic energy is transfered to 
contact energy, it cannot be dissipated anymore, so that 
energy dissipation in the IHS model has to be reduced in 
the presence of multi-particle contacts. This qualitative 
reduction of energy dissipation in dense systems has been 
observed in soft-particle molecular dynamics simulations 

Recently, the transition of a granular gas to its solid 
counterpart has been investigated A general 

model was defined, in which the coefficient of the normal 
restitution is given by 
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if x > x c 
if x < x c 



(1) 



where x is a variable and x c is the respective cut-off value. 
Above the threshold, one has the usual inelastic hard 
sphere mode with constant restitution coefficient, below 
the cut-off an elastic system with r — 1 is assumed. The 
variables proposed were either the time between collisions 
(TC model), the distance travelled since the last collision 
(LC model), or the relative velocity of two particles prior 
to a collision (VC model) . In order to keep the following 
analysis simple, we focus on the special case of a piece- 
wise constant restitution and disregard any continuous 
dependency of r on x |27[] . 

Variants of the general model have been used Jl4|, [H], 
Bjj , \i2\ B3| , mainly to avoid the "inelastic col- 
lapse" , an artefact of the rigid particle model, which al- 
lows an infinite number of collisions to occur within a 
finite time. In the real system this can never occur due 
to the fact that the contacts are finite. The model in Eq. 
(|l|), avoids the inelastic collapse, since the dense parts of 
the system, where the collapse tends to occur first, are 
transformed into elastic regions where the inelastic col- 
lapse is unlikely. Thus the inelastic collapse is replaced 
by static, dense regions of the material. 

The different variants of the cut-off model will be dis- 
cussed separately in the following, because they lead to 
different forms of the collision integral. The general form 
of the energy balance equation is 

^-E = -2I(E,x c ) , (2) 



with the dimensionless time r = (2 / '3) At / 't e{0) , scaled 
by A = (1 — r 2 )/4, and the initial collision rate £_e(0). In 
these units, the energy dissipation rate / is a function of 
the dimensionless energy E = K/K(0) with the kinetic 
energy K, and the cut-off parameter x c . In this represen- 
tation, the restitution coefficient is hidden in the rescaled 
time via A, so that IHS simulations with different r scale 
on the master-curve in the following plots. In the fol- 
lowing, we will extract the classical dissipation rate E 3 / 2 
lUll from /, so that 



I(E,x c ) = J{E,x c )E 3 ' 2 



(3) 



with the correction- function J — > 1 for x c — > 0. Our 
theoretical results will be compared with numerical sim- 
ulations and with previous results Jl6| . For the derivation 
of the dimensionless equation (J2J) from the kinetic theory 
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in its dimensional form, see Refs. |1( 

For the classical IHS model in the HCS, Eq. (|2|) is 
solved by E T = (1 + t)~ 2 , a master curve, independent 
of the coefficient of restitution r and all other system 
parameters. We checked via simulations that different 
r values scale on the same master-curve, as long as no 
clustering is obtained. We will proceed to develop our 
theory in the dimensionless variables and will examine 
in detail the deviations from the classical HCS. 

The velocity cut-off (VC) model can be rationalized 
based on the picture of elasto-plastic particles which do 
not suffer inelastic (plastic) deformation if they collide 
below a certain threshold velocity v c . (In static contact, 
the relative velocity vanishes and thus is automatically 
smaller than v c ), see 14 for a recent application. The 
deviation from the classical inelastic hard sphere HCS, 



J(E,v c )= (l+£ 2 )exp(-£ 2 ) 



(4) 



is obtained from the (cumbersome) computation of the 
collision integral pM , with the nondimensional quantity 



Zmv c 
8K{t) 
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V 2 



(5) 



which relates the critical velocity to the actual mean fluc- 
tuation velocity. The dimensionless cut-off velocity is 
V c = v c /2vt(0)- For v c = and thus £ = 0, the classical 
homogeneous cooling state is recovered, i.e. J(E,0) = 1. 

Event driven numerical simulations pi |l6| are com- 
pared to the numerical solution of our theory in Fig. ||. 
We obtain perfect agreement between theory and sim- 
ulations in the examined range of w c -values. The fixed 
cut-off velocity has no effect when the collision velocities 
are very large, vt 3> v c , but strongly reduces dissipa- 
tion when the relative velocity at collision is comparable 
to or smaller than v c . Thus, in the homogeneous cool- 
ing state, there is no effect initially, but the long time 
behavior changes from the classical decay E oc t~ 2 to a 
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FIG. 2: Energy E plotted against time r for a simulation 
with N = 19683 particles, density v = 0.08245, coefficient of 
restitution ro = 0.99 and different V c as given in the inset. 
The symbols are simulations and the lines give the solution 
of Eqs. d||), using Eqs. (||). 



constant, time independent E cx t , with extremely slow 
approach to the constant. 

The free path cut-off (LC) model was used first by Mc- 
Namara and Young to avoid the inelastic collapse, 
and was extended to its simple hydrodynamic analogon 
expressed in levin-, of the density by Kiinienetskv <•! al. 



23 to describe the gas solid transition caused by the 



compression of a granular gas. The physical idea is that 
particles thath are near to each other - their distance is 
below a certain threshold which can be regarded as some 
surface roughness - are supposed to be in contact with 
each other, so that their contact potential energy cannot 
be dissipated. 

The deviation from the classical inelastic hard sphere 
HCS is J(E,X C ) — J(X C ) cx E°, a constant independent 
of the energy and thus independent of time. Thorough 
calculation |f§] yields 



J(A C ) = exp(-fe A ) 



(G) 



with s\ = \ C (N /V)(Aa) 2 g2a{v) = Ae/y^A, and constant 
k 7.37. This result can be understood, since in the 
homogeneous cooling regime, one has constant density 
and thus constant mean free path, so that a free path 
cut-off model has a time independent effect. Due to its 
lack of interesting new phenomena for the HCS, we will 
not discuss the LC model further. 

The TC model was invented in order to model elas- 
tic material properties, like the "detachment" effect [ p"7| , 
in the framework of the 1HS model. In soft assemblies of 
particles this resembles multi-particle contacts and avoids 
the inelastic collapse in dense IHS systems [jl6|, ^4], p5| ; 
the physical idea behind was discussed in the introduc- 
tion. In technical terms, a collision is elastic if any one 



of two colliding particles had a collision within a time t c 
before the actual time. 

The deviation from the classical HCS is, see flq|, 



J(E,t c ) =exp(*(a:)) 



(7) 



with the series expansion 'i'(x) = — 1.268a; + 0.01682a; 2 — 
0.0005783a; 3 + 0(x 4 ) in the collision integral, with x = 
^t c t E \0)VE = ^t c (0)V~E = ^Fr c |f]. This is close 
to the result V^lm = — 2a;/y / 7r, proposed by Luding and 
McNamara, based on probabilistic mean-field arguments 
[|i"6| psf Here, the argument of the exponential is propor- 
tional to the collision rate t^ 1 cx VE, different from the 
other models, so that J cx exp(v^E). 

Simulation results are compared to the theory in Fig. [|. 
The agreement between simulations and theory is almost 
perfect in the examined range of i c -values, only for large 
deviations from the HCS solution and for large t c -values, 
a few percent discrepancy are observed |29] 

The results can be explained as follows. The fixed cut- 
off time has no effect when the time between collisions 
is very large 3> t c , but strongly reduces dissipation 
when the collisions occur with high frequency t^ 1 ~ i" 1 . 
Thus, in the homogeneous cooling state, there is a strong 
effect initially, but the long time behavior tends towards 
the classical decay E cx t~ 2 . 

Additional simulations with a set of system-sizes and 
for different (also very small ) restitution coefficients will 
be discussed elsewhere Note however, that our con- 
clusions are valid for all system sizes examined and for ar- 
bitrary restitution coefficients before the inhomogeneities 
evolve. 

In summary, a general class of cut-off models was pre- 
sented, aiming towards the enhancement of classical ki- 
netic theory with respect to the realistic behavior of dis- 
sipative particles in the presence of multi-particle inter- 
actions |3Ci|j . Analytical expressions for the collisional 
cooling rate in the energy balance equation of the hy- 
drodynamic equation is provided for the multi-particle 
contacts, evading the singularity of the inelastic collapse. 
Our theoretical results were verified by event-driven nu- 
merical simulations of the HCS and perfect agreement 
was obtained. For realistic material behavior combina- 
tions of the models and also refinements may be nec- 
cessary. Our model, however, leads to a correction of 
the energy dissipation term alone, in the framework of a 
hydrodynamic continuum theory. We regard it thus as 
much simpler than the model proposed in Ref. [^6| that 
also takes the finite contact duration into account. 

The TC model, and to some extent also the other 
models, are based on the assumption that the elastic, 
reversible, potential contact energy of real particles can- 
not be dissipated in the same way as the kinetic energy. 
If one has multi-particle contacts in the system, a lot of 
energy is stored in their contacts - and thus cannot be 
dissipated. 
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FIG. 3: Deviation from the HCS, i.e. rescaled energy E/E T , 
plotted against r for simulations with different r c (0) as given 
in the inset, with rg = 0.99, and N = 8000. Symbols are 
simulation results, the dashed line is the first order correction, 
the solid line results from the third order, and the dotted line 
correspond to the results from [|16|. 



Future interesting work could be the extension of our 
cut-off models to more complicated material laws, e.g., 
introducing some velocity dependent restitution coeffi- 
cient r(v) or contact duration t c (v). In the same spirit, 
the cut-off law can be replaced by continuous functions 
instead of step-functions In addition, the present 

theory should be applied to hydrodynamic models of in- 
homogeneous systems, where the cut-off criterion is a 
function of the position, in order to prove its general ap- 
plicability. As another verification, the model could be 
compared to soft-sphere simulations and experiments. 
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This simplification also is reasonable in the spirit and 
the framework of the kinetic theory and numerical event- 
driven simulations, where changes of the particle veloci- 
ties occur as instantaneous, discontinuous events 
$lm thus neglects non-linear terms and underestimates 
the linear part 

The simulation has to be carefully prepared for the TC 
model in order to achieve good agreement. First, the sys- 
tem is relaxed elastically with r = 1 for several hun- 
dred collisions per particle and the last time of colli- 
sion is saved for each particle. Second, the dissipation 
is switched on and the TC model is activated using the 
saved information about previous contacts. If this infor- 
mation is not used, one observes an artificial initial decay 
of energy in the simulation. 

Only the TC model was used for the detailed discussion 
However, since experimental data are missing, we prefer 
the simple event-based model which is consistent with 
the kinetic theory 



